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Abstract 

Recent  advances  in  computational  materials  science  present  novel  opportunities  for  structure  discovery  and  optimiza¬ 
tion,  including  uncovering  of  unsuspected  compounds  and  metastable  structures,  electronic  structure,  surface  and 
nano-particle  properties.  The  practical  realization  of  these  opportunities  requires  systematic  generation  and  classifi¬ 
cation  of  the  relevant  computational  data  by  high-throughput  methods.  In  this  paper  we  present  Aflow  (Automatic 
Flow),  a  software  framework  for  high-throughput  calculation  of  crystal  structure  properties  of  alloys,  intermetallics 
and  inorganic  compounds.  The  Aflow  software  is  available  for  the  scientific  community  on  the  website  of  the  ma¬ 
terials  research  consortium,  aflowlib.org.  Its  geometric  and  electronic  structure  analysis  and  manipulation  tools  are 
additionally  available  for  online  operation  at  the  same  website.  The  combination  of  automatic  methods  and  user  online 
interfaces  provide  a  powerful  tool  for  efficient  quantum  computational  materials  discovery  and  characterization. 
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1.  Introduction 

A  new  class  of  software  tools  to  assess  material  prop¬ 
erties  has  emerged  from  two  parallel  theoretical  ad¬ 
vancements:  quantum-mechanical  computations  based 
on  density  functional  theory  (DFT),  and  informatics 
data  mining  and  evolutionary  structure  screening  strate¬ 
gies.  Joined  together,  these  methods  make  it  possi¬ 
ble  to  efficiently  screen  large  sets  of  material  structures 
with  many  different  combinations  of  elements,  com¬ 
positions  and  geometrical  configurations.  The  prac¬ 
tical  realization  of  such  schemes  necessitates  a  high- 
throughput  (HT)  approach,  which  involves  setting  up 
and  performing  many  ah  initio  calculations  and  then 
organizing  and  analyzing  the  results  with  minimal  in¬ 
tervention  by  the  user.  The  HT  concept  has  already 
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become  an  effective  and  efficient  tool  for  materials 
discovery  [1,  2,  3,  4,  5,  6,  7,  8]  and  development 
[9,  10,  11,  12,  13].  Examples  of  computational  ma¬ 
terials  HT  applications  include  combinatorial  discov¬ 
ery  of  superconductors  [1],  Pareto-optimal  search  for 
alloys  and  catalysts  [14,  15],  data-mining  of  quan¬ 
tum  calculations  applying  principle-component  analy¬ 
sis  to  uncover  new  compounds  [5,  6,  7,  16,  17],  Kohn- 
anomalies  search  in  ternary  lithium-borides  [18,19,  20], 
and  multi-optimization  techniques  used  for  the  study  of 
high-temperature  reactions  in  multicomponent  hydrides 
[21,22,23]. 

In  its  practical  implementation,  the  HT  approach  uses 
some  sort  of  automatic  optimization  technique  to  screen 
a  library  of  candidate  compounds  and  to  direct  further 
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refinements.  The  library  can  be  a  set  of  alloy  pro¬ 
totypes  [24,  7]  or  a  database  of  compounds  such  as 
the  Pauling  File  [25]  or  the  ICSD  Database  [26,  27]. 
Rather  than  calculating  one  target  physical  quantity  over 
a  large  number  of  structures  and  compositions,  a  key 
philosophy  of  the  HT  method  is  to  calculate  a  priori  as 
many  different  quantities  as  it  is  computationally  feasi¬ 
ble.  Properties  and  property  correlations  are  then  ex¬ 
tracted  a  posteriori. 

This  paper  describes  the  HT  framework  Aflow 
which  we  have  been  developing  over  the  last  decade. 
The  code  and  manuals  describing  its  many  oper¬ 
ation  options  are  freely  available  for  download  at 
aflowlib.org/aflow.html.  Some  of  its  capabilities  may 
also  be  operated  interactively  online  at  a  dedicated  web 
page  aflowlib.org/awrapper.html. 

2.  Software  and  Structure  Database 

The  high-throughput  framework  Aflow,  is  an  assem¬ 
blage  of  software  tools  comprising  over  150,000  lines 
of  C++  code.  It  is  written  for  UNIX  systems  (GNU- 
Linux,  Mac  OSX)  with  the  GNU  suite  of  compilers. 
Its  main  features  are  fully  multi-threaded  and  parallel. 
It  is  designed  to  run  on  top  of  any  software  for  struc¬ 
ture  energy  calculation  and  is  currently  optimized  for 
first-principles  calculations  by  Vienna  Ab  Initio  Simu¬ 
lation  Package  (VASP)  [28]  (porting  to  other  DFT  pack¬ 
ages,  such  as  Quantum  Espresso  [29]  is  underway).  The 
code  offers  several  options  for  running  the  DFT  pack¬ 
age,  on  either  a  single  structure  or  on  sets  of  structures, 
by  searching  through  subfolders  containing  aflow. in  in¬ 
put  files  and  running  those  that  have  not  been  calculated 
yet. 

Aflow  makes  it  possible  to  automatically  calculate  a 
suite  of  physical  observables  over  a  specified  class  or  a 
large  database  of  structures  with  little  human  interven¬ 
tion  to  set  up  the  input  files,  run  the  calculations,  and 
collate  and  plot  the  results.  It  can  also  be  used  to  assist 
in  the  setup  of  standardized  calculations  of  individual 
cases.  To  maximize  the  usefulness  of  the  large  amount 
of  information  produced  by  the  HT  approach,  the  data 
must  be  generated  and  represented  in  a  consistent  and 
robust  manner.  This  is  especially  true  when  the  goal  is, 
as  in  many  materials  science  applications,  to  simultane¬ 
ously  optimize  multiple  properties.  For  instance,  in  cat¬ 
alyst  design  [30]  and  superconducting  materials  devel¬ 
opment  [19,  20],  both  thermodynamics  and  electronic 
structure  are  essential  to  the  effectiveness  of  the  mate¬ 
rial.  If  started  as  a  common  Unix  daemon  through  the 
queue  of  a  computer  cluster,  Aflow  can  generate,  run, 


correct  and  converge  many  calculations  per  day,  with 
minimum  human  input. 

The  backbone  of  the  software  is  the  implementa¬ 
tion  of  ground- state  energy  calculations  to  identify  sta¬ 
ble  and  meta- stable  structures.  In  addition,  it  ad¬ 
dresses  phase  diagram  construction,  electronic  band 
structure,  phonon- spectra  and  vibrational  free-energy, 
surface  stability  and  adsorption  of  contaminants  on  sur¬ 
faces,  electron-phonon  coupling  and  superconductivity, 
local  atomic  environment  expansion  [31],  and  cluster 
expansion.  The  code  is  being  developed  continuously 
with  new  applications  integrated  into  it  periodically. 

One  of  the  most  difficult  challenges  in  HT  database 
generation  is  responding  automatically  to  the  failure  of 
a  calculation.  The  most  common  cause  is  insufficient 
hardware  resources.  Precaution  must  be  taken,  for  ex¬ 
ample,  to  estimate  the  memory  requirement  of  the  tasks 
and  group  jobs  so  that  they  do  not  hamper  each  other 
if  running  on  the  same  node.  The  second  most  com¬ 
mon  cause  of  interruption  is  due  to  runtime  errors  of 
the  ab  initio  calculation  itself.  These  errors  include 
inconsistent  reciprocal  and  k  lattice  meshes,  inconsis¬ 
tent  atom  locations,  ill-defined  geometries,  too  few  elec¬ 
tronic  bands  to  include  all  the  electrons  in  the  system, 
internal  inconsistencies  arising  in  the  convergence  pro¬ 
cess,  and  many  others.  Aflow  is  capable  of  detecting 
most  of  these  problems,  altering  the  input  accordingly 
and  running  again.  This  automatic  restarting  of  runs  is 
achieved  by  diagnosing  the  error  message,  correcting 
the  appropriate  parameters,  and  restarting  the  calcula¬ 
tion. 

As  of  this  writing,  the  Aflow  database  includes 
~400  experimental  prototypes  taken  primarily  from  the 
Pauling  File  [25]  with  some  additions  from  the  ICSD 
[26,  27]  and  the  Navy  Crystal  Lattice  database  [32]. 
The  most  current  list  of  these  experimental  prototypes 
can  be  viewed  by  either  using  the  command-line  ver¬ 
sion  of  Aflow  or  the  online  tool  described  in  Sec.  4. 
Occasionally,  fully  relaxed  calculations  by  Aflow  (cell 
volume  and  shape  and  the  basis  atom  coordinates  inside 
the  cell)  turn  up  as  ground  state  structures  that  do  not 
yet  have  known  experimental  prototypes  [7,  33].  When 
identified,  they  are  also  added  to  the  structure  database. 

The  Aflow  database  also  includes  a  few  million  bcc- 
and  fcc-derived  superstructures  (those  containing  up  to 
20  atoms/cell),  and  a  similar  number  of  hep-derived  su¬ 
perstructures  (all  those  containing  up  to  24  atoms/cell) 
enumerated  formally  using  the  algorithm  of  Refs.  [34] 
and  [35]  (a  derivative  superstructure  is  a  structure  whose 
lattice  vectors  are  multiples  of  those  of  a  parent  lattice 
and  whose  atomic  basis  vectors  correspond  to  lattice 
points  of  the  parent  lattice).  In  practice  only  a  few  hun- 
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Ag2BalSn2_ICSD_25332  (BCT2)  Density  of  States  (States/eV) 


Figure  1:  Left:  Electronic  band  structure  of  Ag2BaSn2,  ICSD  structure  #  25332,  calculated  along  the  Brillouin  Zone  path  shown  in  the  inset.  Right: 
The  electronic  density  of  states  of  the  same  structure.  Both  automatically  generated  using  the  Aflow  framework. 


dred  of  these  structures  are  usually  calculated  for  each 
alloy  system,  but  they  are  useful  in  constructing  cluster 
expansions  [36,  37,  38,  39]  which  can  be  employed  in 
synergy  with  the  HT  calculations  to  examine  structural 
configurations  more  thoroughly  than  is  possible  with  a 
purely  HT  approach  [17]. 

3.  Basic  operation 

At  the  beginning  of  a  calculation,  the  starting  struc¬ 
ture,  or  structures,  are  selected  from  the  database.  The 
database  structural  description  contains  the  lattice  pa¬ 
rameters  and  atomic  positions  of  the  prototype  com¬ 
pound.  Aflow  then  adjusts  the  lattice  parameters  (us¬ 
ing  a  weighted  average  of  the  pure  element  volumes) 
for  the  specified  elements  and  creates  an  input  file  with 
all  the  necessary  parameters  for  relaxations,  static  and 
band  structure  runs.  This  entire  procedure  requires  only 
the  label  of  a  structure  in  the  Aflow  database.  It  may 
be  performed  with  equal  ease  for  only  one  structure 
with  a  specific  composition  or  for  a  large  set  of  struc¬ 
tures  with  many  combinations  of  elements,  with  a  sin¬ 
gle  command.  For  multiple  structures  or  compositions 
this  command  creates  a  set  of  subfolders,  each  includ¬ 
ing  a  single  input  file,  through  which  Aflow  runs  auto¬ 
matically  until  all  structures  and  compositions  are  cal¬ 
culated. 


A  HT  computational  framework  must  contain  a  gen¬ 
eral,  reliable,  and  standardized  electronic  structure  anal¬ 
ysis  feature.  For  example,  it  must  automatically  deter¬ 
mine  the  Brillouin  Zone  (BZ)  integration  path  for  the 
14  Bravais  Lattices  (with  their  24  Brillouin  zones)  [40] 
and  change  the  basis  (lattice  vectors)  into  a  standard¬ 
ized  basis  so  that  data  can  be  compared  consistently  be¬ 
tween  different  projects.  Although  BZ  integration  paths 
have  appeared  in  the  literature  in  the  last  few  decades 
[41,  42,  43,  44,  45],  a  standardized  definition  of  the 
paths  for  all  the  different  cases  has  been,  to  the  best 
of  our  knowledge,  missing,  and  mistakes  in  the  litera¬ 
ture  for  less-common  Brillouin  zones  are  not  uncom¬ 
mon.  This  component  of  Aflow  was  discussed  in  detail 
in  a  separate  publication  [40]. 

Aflow  computes  structure  total  energies  and  elec¬ 
tronic  band  structures  using  vasp  with  pseudotentials 
and  accuracy  chosen  by  the  user.  As  a  default,  Aflow 
employs  projector  augmented  wave  (PAW)  pseudopo¬ 
tentials  with  GGA  exchange  correlation  functionals  [46, 
47]  as  parameterized  by  Perdew-Burke-Ernzerhof  [48]. 
Each  structure  is  fully  relaxed  twice  with  a  convergence 
tolerance  of  1  meV/atom  using  dense  grids  of  3,000- 
6,000  k-points  per  reciprocal  atom.  At  the  beginning  of 
these  structural  relaxation  steps,  a  spin-polarized  calcu¬ 
lation  is  performed  for  all  structures.  If  the  magnetiza¬ 
tion  is  smaller  than  0.025  yu^/atom,  the  spin  is  turned  off 
for  the  next  relaxations  and  subsequent  calculations  to 
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enhance  the  calculation  speed.  After  this,  the  structure 
is  changed  into  a  standard  form  [40],  and  another  elec¬ 
tronic  relaxation  is  performed  using  fixed  coordinates 
for  the  lattice  vectors  and  atomic  positions.  This  static 
run  is  implemented  with  a  much  denser  grid  of  10,000- 
30,000  points  to  get  accurate  charge  densities  and  den¬ 
sity  of  states  for  the  calculation  of  the  band  structure. 
The  Monkhorst-Pack  scheme  [49]  is  employed  in  the 
grid  generation  except  for  hexagonal,  fee  and  rhombo- 
hedral  systems  in  which  T-centered  grid  is  selected  for 
faster  convergence.  This  process  is  performed  in  one 
calculation,  with  the  single  input  file  Aflow  creates  as 
described  above.  At  the  completion  of  such  a  calcu¬ 
lation  Aflow  invokes  appropriate  Matlab  or  Gnuplot 
scripts  for  data  analysis  and  visualization.  All  these 
steps  are  usually  performed  automatically,  but  the  user 
also  has  the  option  to  operate  them  through  the  web 
interface.  An  example  of  a  computed  electronic  band 
structure  and  standard  BZ  path  appears  in  Figure  1 . 

Aflow ’s  capability  to  continuously  search  subfolders 
for  calculations  to  run  is  not  limited  to  DFT  calcula¬ 
tions.  An  “alien”  mode  is  also  implemented,  which 
allows  Aflow  to  execute  other  tasks  in  a  high  through¬ 
put  fashion.  For  instance,  the  many  thousands  of  grand 
canonical  Monte  Carlo  calculations  used  in  a  recent 
surface  science  absorption  project  [50,  51]  were  di¬ 
rected  and  performed  by  Aflow.  In  addition,  Aflow  is 
equipped  with  options  to  run  commands  or  scripts  be¬ 
fore  and  after  the  main  program  is  performed  in  each 
folder.  This  allows  Aflow  to  generate  input  files  on 
the  fly  depending  on  the  results  of  different  calculations, 
so  that  ad-hoc  optimization  can  be  implemented  by  the 
user.  It  also  improves  the  flexibility  of  recovery  from 
a  crash  or  an  unconverged  run  and  increases  the  overall 
versatility  and  throughput  of  the  calculation. 

4.  Structure  analysis  tools 

Aflow  offers  structure  analysis  and  manipulation 
tools  which  are  also  useful  to  users  who  do  not  need 
to  perform  HT  calculations  or  to  create  databases 
for  datamining.  These  users  may  prepare  standard 
unit  cell  input  files  and  extract  the  appropriate  k- 
points  path  using  the  command  version  of  Aflow, 
called  Aconvasp,  or  the  online  interface  on  our  website 
aflowlib.org/awrapper.html,  (shown  in  Figure  2).  Runs 
should  then  be  performed  according  to  the  following 
protocol:  Unit  cells  must  first  be  reduced  to  standard 
primitive,  then  appropriately  relaxed  by  Aflow.  Before 
the  static  run,  the  cells  should  be  reduced  again  to  stan¬ 
dard  primitive,  since  symmetry  and  orientation  might 
have  changed  during  the  relaxation.  The  user  should 


then  perform  the  static  run  and  then  project  the  eigen¬ 
values  along  the  directions  which  are  specified  in  the  “k- 
path”  option.  If  the  user  is  running  Aflow  and  vasp,  the 
web  interface  can  also  prepare  the  input  file,  aflow. in, 
which  will  perform  all  the  mentioned  tasks.  The  conver¬ 
sion  and  analysis  operations  may  be  carried  on  a  struc¬ 
ture  file  stored  in  the  database,  or  supplied  by  the  user 
in  the  input  box  Input  POSCAR.  The  operations  imple¬ 
mented  on  the  web  interface  are  the  following: 

•  Normal  Primitive;  Generates  a  primitive  cell  in  the 
most  compact  form  (not  necessarily  unique). 

•  Standard  Primitive;  Chooses  a  primitive  cell  such 
that  the  Wigner-Seitz  cell  defined  by  the  recipro¬ 
cal  vectors  coincides  with  one  of  the  possible  24 
Brillouin  zones  [42,  41,  40]. 

•  Standard  Conventional;  Generates  a  conventional 
unit  cell  (not  necessarily  primitive). 

•  Minkowski  lattice  reduction;  Generates  a 
maximally-compact  cell  (not  necessarily  unique). 

•  Niggli  Standardized  form;  Generates  a  cell  that 
conforms  to  the  Niggli  standard  form  [52]. 

•  WYCKOFF-CAR/ABCCAR  to  POSCAR;  Gener¬ 
ates  a  POSCAR  file  from  a  structure  file  containing 
the  standard  crystallographic  information. 

•  POSCAR  to  ABCCAR;  Generates  a  structure  file 
containing  lattice  parameters  (rather  than  lattice 
vectors)  from  a  POSCAR  file. 

•  Bring  atoms  in  the  cell;  Remaps  atomic  coordi¬ 
nates  to  lie  inside  the  unit  cell. 

•  Cartesian  coordinates;  Converts  the  atomic  coor¬ 
dinates  in  a  structure  file  to  Cartesian  coordinates 
(from  lattice  coordinates). 

•  Fractional  coordinates;  Converts  atomic  positions 
from  fractional  coordinates  (i.e,  lattice  or  direct  co¬ 
ordinates)  to  Cartesian  coordinates. 

•  Data  and  Extended  Data;  Generates  lattice  param¬ 
eters  (cell  lengths  and  angles)  and  other  informa¬ 
tion  (volume,  reciprocal  lattice,  symmetry  infor¬ 
mation,  etc.). 

•  Symmetry  Information;  Lattice  Type  of  the  crystal 
and  the  lattice,  Pearson  symbol,  Space  group  and 
Wyckoff  positions,  Point  Group  Lattice  Matrices, 
Point  Group  Crystal  Matrices  and  factor  Group 
Crystal  Matrices  and  Translations. 
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O  aflow/aconvasp  VERSION  O  aflow  HELP  O  aconvasp  HELP  O  apennsy  HELP 
O  Normal  Primitive  (fast) 

O  Standard  Primitive  (slow)  [DOI:  10.1016/j.commatsci.2010.05.010] 

O  Standard  Conventional  (slow)  [DOI:  10.1016/j.commatsci.2010.05.010] 

O  Minkowski  lattice  reduction  O  Niggli  Standardized  form 

O  WYCKOFF-CAR/ABCCAR  to  POSCAR  (read  aconvasp.pdf  for  definition  of  wyckcar/abccar) 

O  POSCAR  to  ABCCAR 

O  Bring  atoms  in  the  cell  O  Cartesian  coordinates  O  Fractional  coordinates 

O  Data  (a, b,c, alpha, beta, gamma)  O  Extended  Data  (a, b,c, alpha, beta, gamma  lattices,  point  group) 

O  Lattice  Type  of  the  crystal  O  Lattice  Type  of  the  lattice  O  Pearson  symbol 
O  Space  Group  (PLATON)  ©  New  Space  Group  (FINDSYM)  O  New  Findsym  output 
O  Point  Group  Lattice  Matrices  O  Point  Group  Crystal  Matrices 

O  Factor  Group  Crystal  Matrices  and  Translations  O  Identical  atoms  O  Sites  Simmetries 
O  Kpath  in  the  reciprocal  space  for  band  structure  calculations  [DOI:  10.1016/j.commatsci.2010.05.010] 
(with  respect  to  the  reciprocal  vectors  generated  from  the  standard  primitive  lattice  vectors). 

O  New  Interstitial  Positions  (online  up  to  8  atoms/cell)  [DOI:  10. 1103/PhysRevB. 79. 134203] 

O  List  of  prototypes  of  AFLOW 
O  List  of  Bands  (reference  needed) 

O  List  of  Libraries  (reference  needed) 

O  Test  for  development 

reference  (if  required,  otherwise  leave  it  empty) 


(be  patient) 
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Figure  2:  Input  interface  for  the  structure  analysis  tools  of  Aflow. 
Users  may  apply  these  tools  to  a  structure  included  in  the  database  or 
enter  their  own  structure  in  the  input  box.  The  tools  include  symme¬ 
try  analysis,  format  conversions,  and  transforming  structures  between 
equivalent  representations  (primitive  versus  conventional,  etc.). 


•  Identical  atoms  and  site  symmetries;  Identifies 
which  atoms  in  a  structure  correspond  to  the  same 
WyckofF  positions. 

•  K-path  in  the  reciprocal  space;  Provides  direc¬ 
tions  in  k-space  needed  for  band  structure  calcu¬ 
lations  and  provides  a  picture  of  the  corresponding 
Brillouin  zone  and  path. 

•  Interstitial  Cages  Positions;  provides  all  the  geo¬ 
metric  locations  of  interstitials  in  the  structure. 

The  user  interface  also  provides  access  to  the  soft¬ 
ware  manuals,  under  the  buttons  marked  aflow  HELP, 
aconvasp  HELP  and  apennsy  HELP,  and  to  the  current 
list  of  experimental  prototypes  in  the  Aflow  database 
under  List  of  prototypes  of  AFLOW. 

5.  Example  Applications 

5.1.  Vibration  spectra  and  free  energy 

Aflow  can  calculate  phonon  dispersion  curves  us¬ 
ing  three  different  approaches:  the  direct  force  constant 
method,  following  the  definition  of  Ref.  [54,  55]  and 
description  of  Ref.  [56],  the  linear  response  method 


Figure  4:  The  supercells  of  (001),  (110),  and  (111)  slabs  of  fee  binary 
Llo  structure  built  by  Aconvasp. 


for  PAWs  [57],  and  the  frozen  phonon  method  as  imple¬ 
mented  in  the  Frozsl  code  of  H.  Stokes  and  L.  Boyer 
[58,  59,  60]. 

In  the  first  method,  the  atomic  displacements  neces¬ 
sary  to  fully  determine  the  dynamical  matrix  are  pre¬ 
scribed  automatically,  the  forces  are  calculated  by  vasp 
and  mapped  onto  symmetry-equivalent  directions  by 
Aflow.  The  method  is  general  and  can  treat  low- 
symmetry  structures  as  easily  as  high-symmetry  sys¬ 
tems,  minimizing  user- time  overhead.  To  reduce  com¬ 
putational  efforts,  the  symmetry  of  the  supercells  is  fully 
employed  and  only  the  unique  atoms  are  distorted  along 
the  minimum  independent  directions  required  for  the 
construction  of  the  force  constant  matrix.  When  neces¬ 
sary,  we  use  linear  response  [57]  for  the  non-analytical 
part  of  the  dynamical  matrix  [61,  62],  to  reproduce  the 
correct  LO-TO  splitting  (splitting  between  longitudinal 
and  transverse  optical  phonon  frequencies).  For  the  lin¬ 
ear  response  method  as  implemented  in  vasp,  Aflow 
sets  up  the  input  files,  runs  and  converges  the  calcula¬ 
tion  automatically.  Aflow  uses  the  generated  disper¬ 
sion  curves  to  calculate  the  vibrational  contribution  to 
the  free  energy,  entropy  and  specific  heat  (see  example 
in  Figure  3).  The  whole  task  is  performed  by  the  APL 
library  of  Aflow. 

In  the  frozen-phonon  approach,  Aflow-Frozsl,  the 
user  sets  up  an  aflow. in  input  file  containing  instruc- 
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Figure  3:  Top:  Phonon  dis¬ 
persion  curves  for  face-centered 
cubic  rhodium.  Solid  lines 
are  computed  frequencies,  dots 
are  experimental  results  [53]. 
The  phonon  density  of  states 
is  shown  on  the  right  hand 
side.  Bottom:  Vibrational 
free  energy,  vibrational  entropy, 
and  specific  heat  calculated  au¬ 
tomatically  using  the  Aflow 
framework. 


tions  in  the  Frozsl  format  [60].  To  perform  the  cal¬ 
culation  Aflow  repeatedly  calls  Frozsl  to  generate  all 
the  irreducible  representations  of  the  phonons,  calculate 
the  various  geometries  and  the  phonon  spectrum  in  the 
requested  k  points  or  paths.  The  code  is  multi- threaded: 
the  various  irreducible  representations  can  be  tackled 
simultaneously  in  a  computer  cluster  while  an  Aflow 
daemon  waits  for  their  completion  to  produce  the  spec¬ 
trum.  For  instance  the  several  thousand  calculations  of 
Ref.  [63]  were  rapidly  addressed  by  this  method. 

We  chose  to  implement  all  three  approaches  for  con¬ 
venience.  While  the  first  gives  more  information  in  the 
whole  Brillouin  Zone,  the  second  and  the  third  can  be 
used  to  rapidly  scan  spectra  at  particular  k  points,  such 
as  when  searching  for  Kohn  anomalies  [20] .  Currently, 
the  Frozsl/Findsym  source  packages  are  maintained  by 
the  Aflow  team.  They  are  included  in  the  Aflow  dis¬ 
tribution  and  have  been  modified  to  compile  with  the 
GNU  suite. 

5.2.  Design  of  high-index  surfaces  in  complex  multi- 
component  compounds 

Surface  segregation,  adsorption,  chemisorption,  cat¬ 
alytic  reaction  and  other  surface  phenomena  are  of  ulti¬ 
mate  importance  in  many  areas  of  modern  technology, 
including  catalysis,  corrosion  protection,  batteries,  elec¬ 
tronics  etc.  In  most  theoretical  studies  of  surface  effects, 
it  is  necessary  to  know  the  exact  positions  of  surface  and 
subsurface  atoms.  In  cases  of  low-index  surfaces  and 
simple  compounds,  the  construction  of  a  surface  can  be 


done  manually.  However,  the  technologically  interest¬ 
ing  high-index  surfaces  [64,  65]  and  multicomponent 
compounds  with  complex  underlying  crystal  lattices  re¬ 
quire  an  automated  tool. 

The  input  file  used  for  the  surface  construction  de¬ 
scribes  the  initial  bulk  crystal  structure,  denoted  by  the 
corresponding  Aflow  database  label.  The  structure  may 
have  an  arbitrary  complex  crystal  lattice  with  any  num¬ 
ber  of  atom  types.  The  generated  surface  file  con¬ 
tains  the  complete  data  (both  direct  and  Cartesian  co¬ 
ordinates)  for  the  constructed  surface  supercell.  The 
Miller  indices  (hkl)  are  based  on  the  Bravais  lattice  that 
corresponds  to  the  unit  cell  of  the  bulk  crystal  struc¬ 
ture.  An  example  of  (001),  (110),  and  (111)  super¬ 
cell  slabs  of  the  fee  binary  Llo  structure  are  shown  in 
Figure  4.  Aflow  can  easily  build  a  supercell  contain¬ 
ing  a  predetermined  number  Nf  of  atomic  planes  with 
given  Miller  indices  (hkl)  of  a  given  crystal  structure. 
A  supercell  may  also  contain  a  number  Ne  of  “empty” 
(MZ)-planes,  without  actual  atoms.  In  the  last  case,  a 
periodically  repeated  supercell  represents  the  slabs  sep¬ 
arated  by  empty  space.  Such  slab  construction  allows 
one  to  study  surface  effects,  e.g.  surface  energy,  sur¬ 
face  stress,  and  surface  segregation,  from  first  princi¬ 
ples  [66,  67,  68,  69,  70,  65].  The  obtained  surface  data 
may  also  be  effectively  used  in  simulation  of  nanopar¬ 
ticles  [70,  65].  The  variation  of  Nf  and  Ne  is  im¬ 
portant  for  studying  convergence  with  respect  to  slab 
width  and  inter  slab  distance,  respectively.  Within  (hkl)- 
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planes,  the  constructed  supercell  consists  of  p\  x/?2  two- 
dimensional  primitive  unit  cells,  where  pu  P2  =  1,2,... 
are  specified.  Large  p\  and  P2  are  used  to  introduce  sur¬ 
face  defects,  e.g.  vacancy  or  solute  atoms,  far  enough 
from  each  other  to  diminish  their  interaction.  For  in¬ 
stance,  all  the  surfaces  studied  in  Ref.  [65]  were  pro¬ 
duced  with  this  method. 

5.3.  Nanoparticle  generation 

Starting  from  an  input  crystal  structure,  a  specified 
radius  and  required  separation,  Aconvasp  generates  a 
structure  file  of  a  nano-particle  of  that  radius  made  of 
the  same  lattice  and  separated  as  specified  from  its  near¬ 
est  neighbors.  The  origin  of  the  particle  may  be  set  to 
the  origin  of  the  unit  cell,  an  atom  in  the  crystal  struc¬ 
ture  or  any  point  in  Cartesian  or  fractional  coordinates. 
The  radius  and  separation  distance  are  in  A.  This  option 
replaces  the  cumbersome  manual  generation  of  nano¬ 
particle  structure  files,  usually  employed  in  studies  of 
such  particles,  and  should  enable  investigation  of  large 
sets  of  nano-particles  in  a  high-throughput  fashion 


Figure  5:  Nano-particles  based  on  the  LI3  structure  built  by  Aconvasp, 
with  radii  of  5A  (left),  8A  (center)  and  12A  (right). 


located  in  any  of  those  close  positions  could  deform  the 
nearby  local  atomic  environment  and  relax  to  the  same 
final  location.  Hence,  the  number  of  symmetrically  in¬ 
equivalent  cages  can  be  further  reduced  by  agglomera¬ 
tion  of  such  a  set  of  positions  into  a  single  interstitial 
site,  upon  insertion  of  interstitial  atoms. 

Given  the  host  cell  geometry  and  the  interstitial 
species,  the  occupation  of  the  final  irreducible  cages  can 
be  automatically  simulated  by  Aflow,  which  calculates 
their  energies,  entropy  and  solubility  in  the  small  con¬ 
centration  limit  [71,  72].  These  operations  are  imple¬ 
mented  in  a  multi-threaded  manner  to  accelerate  the  cal¬ 
culation  in  a  multi-core  environment.  Figure  7  presents 
an  example  of  the  output  of  a  multi-threaded  interstitial 
position  search  performed  online  through  our  servers 
(aflowlib.org/awrapper.html),  for  the  AgZr3-Ll3  struc¬ 
ture,  using  the  interface  of  Figure  2.  Here  Aflow  finds 
14  cages:  6  octahedral  and  8  tetrahedral.  Reduction 
through  the  space  group  leads  to  only  5  irreducible 
positions,  two  octahedral  positions  with  r  =  2.2225 
A(AgiZr5  and  Ag2Zr4)  and  three  tetrahedral  positions 
with  r  =  1.9248  A(AgiZr3,  Ag2Zr2  and  Zr4).  The 
equivalence  and  the  multiplicity  per  cell  are  provided  in 
the  text.  Figure  6  shows  these  octahedral  and  tetrahedral 
interstitials  in  the  standard  conventional  representation 
of  this  structure  Ll3.  For  clarity,  the  conventional  0P8 
structure  is  shown,  instead  of  the  primitive  0S8.  The  two 
unique  octahedrons  are  blue  and  the  three  unique  tetra- 
hedra  are  in  pink.  The  interstitial  atom  is  in  the  center 
of  the  poly  topes. 


5.4.  Topological  identification  of  interstitial  sites 

An  additional  feature  of  Aflow  allows  users  to  iden¬ 
tify  interstitial  sites  inside  any  crystallographic  struc¬ 
ture.  Using  as  input  a  file  specifying  the  atomic  posi¬ 
tions  in  the  structure  unit  cell,  it  returns  the  location,  ra¬ 
dius  and  coordination  of  each  site.  The  algorithm  iden¬ 
tifies  quadruplets  of  non-coplanar  atoms,  where  the  first 
atom  belongs  to  the  unit  cell  and  the  others  are  sepa¬ 
rated  by  less  than  the  longest  diagonal  of  the  unit  cell. 
Cages  are  defined  by  spheres  touching  all  four  atoms  of 
the  quadruplet  that  do  not  contain  any  atoms.  An  inter¬ 
stitial  position  is  found  if  a  center  of  a  cage  is  inside  the 
unit  cell.  By  considering  all  the  possible  combinations, 
symmetrically  inequivalent  interstitials  can  be  identi¬ 
fied  through  calculation  of  their  site  symmetry  (with  the 
space  group  of  the  unit  cell).  Note  that  in  unit  cells  with 
complex  arrangements,  many  of  the  interstitial  positions 
can  be  extremely  close.  Thus,  adjacent  interstitial  atoms 
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Figure  6:  6-coordinated  (blue)  and  4-coordinated  (pink)  interstitial 
sites  in  the  LI3  structure,  found  in  a  topological  search  by  Aconvasp. 
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S  CALCULATION  -  S 


REFERENCE  STRUCTURE  BELOW 


AgZr/11  -  (11)  -  FCC  [AB3 ]  (11)  (htqc  library) 


UNIQUE=14  -  cage4  radius  max=2 . 2225985  -  cage4  radius  min=l . 92482 68  - 

REDUCIBLE:  14 


Direct (5)  [A1B1C3] 

0.000000  0.000000  0.000000  Ag 

0.000000  0.625000  0.875000  H 

0.500000  0.000000  0.500000  Zr 

0.500000  0.500000  0.000000  Zr 

0.000000  0.500  0  00  0.500000  Zr 

[VASP_POSCAR_MODE_EXPLICIT] STOP . cages_T3 

[VASP_POSCAR_MODE_EXPLICIT] START . cages_T4 


[VASP_POSCAR_MODE_EXPLICIT] START . cages_T5 

AgZr/11  -  (11)  -  FCC  [AB3]  (11)  (htqc  library)  R=1 . 92483  T=5 


[Zr (1) ,Ag(0) ,Zr (2) ,Ag(0) ] 


Figure  7:  Output  of  a  multi-threaded  search  of  interstitial  cages 
in  the  LI3  structure,  showing  two  octahedral  positions  with  r  - 
2.2225 A(AgiZr5  and  Ag2Zr4)  and  three  tetrahedral  positions  with 
r  =  1.9248A(AgiZr3,  Ag2^2  and  Zr4).  The  cell  geometries  includ¬ 
ing  these  interstitials  are  given  in  the  bottom  part  of  the  printout. 


5.5.  Apennsy:  automatic  analysis 

The  results  of  high-throughput  Aflow  runs  are  an¬ 
alyzed  and  manipulated  for  further  processing  by  the 


Apennsy  module.  The  basic  operation  of  this  module 
is  the  production  of  a  system-specific  data  file  that  in¬ 
cludes  information  on  all  the  computed  structures,  their 
energies,  the  input  and  output  (after  relaxation)  crys¬ 
tal  prototype  and  space  group,  and  the  binary  convex 
hull.  It  requires  as  input  the  output  files  generated  by 
the  Aflow  run  for  the  computed  system. 

For  example,  in  a  binary  alloy  system,  the  convex  hull 
at  zero-temperature  is  used  to  create  the  phase  diagram. 
To  construct  the  convex  hull,  Apennsy  finds  the  mini¬ 
mum  energy  structure  at  each  composition,  and  calcu¬ 
lates  their  formation  enthalpies 


Hf  =  E(x)-xEA-(l-x)EB,  (1) 

where  Ea  and  EB  are  the  energies  of  the  pure  elements 
and  E(x)  is  the  energy  of  the  alloy  with  concentration 
v  of  element  A.  Apennsy  then  identifies  the  extremal 
structures,  i.e.,  those  that  lie  below  the  tie-lines  connect¬ 
ing  their  neighbors  at  adjacent  compositions.  These  ex¬ 
tremal  structures  are  the  stable  structures  in  the  binary 
system  and  constitute  its  convex  hull.  Apennsy  plots  the 
convex  hull  automatically  using  Gnuplot  or  Matlab. 

The  module  can  also  be  used  to  print  the  computed 
data  in  formats  appropriate  for  use  as  input  for  other 
codes.  For  example,  it  can  print  all  the  bcc  (or  fee  or 
hcp)-based  structures  with  the  total  energy  per  unit  cell 
or  energy  per  atom.  This  is  useful  for  cluster  expan¬ 
sion  codes  such  as  uncle  [39].  Another  format  is  used 
for  automatic  miscibility  determination.  All  these  data 
files  can  be  generated  off-line  by  operating  Apennsy  on 
a  library  of  input  files,  or  on-line  by  selecting  a  spe¬ 
cific  system  and  the  requested  option  at  the  Aflow  web¬ 
site  aflowlib.org  (“binary  alloy  library”  entry).  Addi¬ 
tional  options  of  Apennsy  provide  data  that  is  external 
to  the  structure  calculations,  such  as  predictions  of  sys¬ 
tem  miscibility  following  Miedema  and  Hume-Rothery 
mixing  rules  [73,  74]. 


6.  Future  developments:  High-throughput  hybrid 

functionals  calculations 

The  local  (LDA)  and  semi-local  (GGA)  approxima¬ 
tions  to  the  DFT,  expanded  with  “+U”  techniques  [75], 
have  been  incredibly  successful  in  describing  a  wide  va¬ 
riety  of  materials  properties,  particularly  those  related 
to  the  ground  state.  In  some  cases,  however,  the  draw¬ 
backs  of  these  approximations  present  an  insurmount¬ 
able  obstacle.  One  of  the  most  illustrative  areas  is 
the  energetics  of  defects  in  wide  band-gap  semiconduc¬ 
tors.  For  example,  in  the  case  of  oxygen  vacancies  in 
ZnO,  the  defect  formation  energies  and  thermodynamic 


transition  levels  obtained  by  different  LDA/GGA  ap¬ 
proaches  vary  by  several  eV  [76].  Similar  problems  ex¬ 
ist  in  DFT  calculations  of  optical  properties  of  insulators 
and  semiconductors,  properties  of  correlated  materials, 
reaction  energetics  and  adsorption  of  small  molecules 
[77,  78].  An  affordable  and  accurate  approach  to  over¬ 
come  these  problems  does  not  exist  yet.  The  state  of 
the  art  many-body  methods,  such  as  self-consistent  GW 
[79]  exhibit  excellent  agreement  with  experiments,  but 
suffer  from  high  computational  cost.  The  hybrid  func¬ 
tional  methods  [80]  represent  a  pragmatic  compromise 
between  the  local  and  nonlocal  many-body  approaches. 
The  results  of  hybrid  functional  calculations  are  often 
in  excellent  agreement  with  computationally  expensive 
cutting  edge  many-body  methods,  yet  favorable  scaling 
allows  routine  calculations  of  a  few  hundred  atoms  in 
the  unit  cell  [81]  . 

In  a  hybrid  functional  calculation  the  LDA  exchange 
correlation  part  of  the  density  functional  is  mixed  with 
a  Fock-type  exchange  part  in  varying  proportions.  For 
example,  the  popular  Perdew-Burke-Ernzerhof  hybrid 
functional  (PBEO)  contains  25%  of  the  exact  exchange, 
75%  of  GGA  exchange,  and  100%  of  GGA  correla¬ 
tion  energy  [82].  In  a  related  Heyd-Scuseria-Emzerhof 
(HSE06)  functional  the  Fock  exchange  interactions  are 
separated  into  long-  and  short-range  parts.  The  short 
range  part  includes  25%  of  exact  exchange  and  75%  of 
semi-local  GGA  exchange,  while  the  long-ranged  part  is 
replaced  with  an  approximate  semi-local  GGA  expres¬ 
sion  [83].  The  splitting  is  accomplished  by  introducing 
the  screening  of  the  exchange  interactions  (similar  to  the 
screened  exchange  approach)  with  optimal  screening 
length  of  approximately  7-10  A.  This  approach  presents 
a  middle  case  between  the  GGA  calculation  (no  exact 
exchange)  and  the  PBEO  hybrid  functional  (all  exchange 
is  long  ranged).  The  HSE06  functional  eliminates  some 
unphysical  features  of  the  exact  exchange  approach.  It 
also  has  the  computational  advantage  of  better  conver¬ 
gence  of  the  long-ranged  exchange  part,  since  screening 
out  the  long  range  exact  exchange  greatly  reduces  the 
cost  of  calculating  the  non-local  exchange  interaction. 
One  also  often  finds  it  useful  to  tune  the  amount  of  exact 
exchange  for  a  particular  material  (the  so  called  cr-tuned 
hybrid  functionals)  in  order  to  obtain  the  best  agreement 
of  computed  band-gap  with  experiment.  This  is  par¬ 
ticularly  useful  in  cases  where  the  correct  value  of  the 
band-gap  is  critical.  For  example,  in  the  calculations  of 
defect  energetics,  tuning  the  band-gap  of  the  host  mate¬ 
rial  in  the  bulk  and  then  running  calculations  for  super¬ 
cells  containing  impurities  and  defects  has  been  shown 
to  produce  results  in  very  good  agreement  with  experi¬ 
ment  [84]  . 


The  hybrid  functional  calculations  are  computation¬ 
ally  more  demanding  than  LDA/GGA  calculations,  due 
to  the  need  to  fully  calculate  the  non-local  Fock  ex¬ 
change  integrals.  For  example,  a  bulk  ZnO  HSE06  cal¬ 
culation  (8x8x6  T-centered  k-point  mesh,  with  40 
irreducible  k-points,  400  eV  cutoff)  using  vasp  takes 
approximately  5  hours  on  8  CPU’s.  The  same  set  up 
for  the  GGA/PBE  (+U)  calculation  takes  less  than  a 
minute.  The  memory  requirements  in  this  case,  for  a 
HSE06  calculation,  are  also  increased  by  a  factor  of  five 
compared  to  those  of  GGA.  The  scaling  of  the  computa¬ 
tional  cost  varies  for  different  hybrid  functionals,  and  is 
approximately  0(N2-5)  for  small  system  sizes,  and  lin¬ 
ear  for  HSE06  beyond  15  A  and  for  PBEO  beyond  100 
A.  Using  plane  waves  the  scaling  of  hybrid  functional 
methods  is  (N^^N^N^Tlog (  Afft ) ,  which  is  approx¬ 
imately  linear  with  the  number  of  atoms  in  the  bulk 
[81].  The  standard  (semi-)  local  DFT(+U)  methods  us¬ 
ing  plane  waves  normally  scale  as  0(N3). 

6. 1 .  Framework  of  DFT +  U  and  Hybrid  Functional  cal¬ 
culations 

HT-DFT+U:  Within  the  high-throughput  frame¬ 
work,  Aflow  currently  employs  the  DFT+U  technique 
and  allows  the  user  to  choose  the  appropriate  parame¬ 
ters.  As  default,  Aflow  takes  advantage  of  the  Dudarev 
[85]  formalism  within  GGA+U.  The  values  of  U  are 
listed  in  Table  1. 


Table  1:  Values  of  C/eff  parameters  in  eV  for  the  Dudarev 
GGA+U  approach  implemented  in  Aflow.  From  Refs.  [86, 
40]  and  references  therein. 


Ti 

V 

Cr 

Mn 

Fe 

Co 

Ni 

Cu 

Zn 

Ga 

4.4 

2.7 

3.5 

4.0 

4.6 

5.0 

5.1 

4.0 

7.5 

3.9 

Nb 

Mo 

Tc 

Ru 

Rh 

Pd 

Cd 

In 

2.1 

2.4 

2.7 

3.0 

3.3 

3.6 

2.1 

1.9 

Ta 

W 

Re 

Os 

Ir 

Pt 

2.0 

2.2 

2.4 

2.6 

2.8 

3.0 

La 

Ce 

Pr 

Gd 

Nd 

Sm 

Eu 

Tm 

Yb 

Lu 

7.5 

6.3 

5.5 

6.0 

6.2 

6.4 

5.4 

6.0 

6.3 

3.8 

HT-DFT+Hybrid:  We  are  currently  expanding  the 
Aflow  framework  to  perform  HSE06  calculations,  and 
we  present  some  preliminary  results.  The  extension  is 
not  trivial:  the  calculations  of  electronic  structure  using 
hybrid  functionals  differ  from  the  standard  LDA/GGA 
ones,  with  or  without  “+U”.  One  cannot  compute  the 
bands  in  the  familiar  non-selfconsistent  way  of  the 
LDA/GGA  because  non-local  exchange  is  not  deter¬ 
mined  by  the  pre-computed  charge  density.  Therefore  to 
obtain  the  eigenvalues  for  strings  of  k-points  the  follow¬ 
ing  recipe  is  being  implemented.  First,  Aflow  performs 
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a  standard  LDA/GGA  calculation.  Second,  a  hybrid 
functional  run  is  performed  starting  from  the  LDA/GGA 
wavefunctions  on  the  same  k-point  mesh  and  energy 
cutoff  chosen  by  Aflow.  The  number  of  bands  is  dy¬ 
namically  adjusted  to  achieve  full  convergence.  Then, 
to  facilitate  efficient  use  of  computational  resources,  it 
is  further  adjusted  to  include  only  a  few  bands  over  the 
highest  occupied  state.  Third,  the  electronic  structure 
is  computed  by  performing  a  hybrid  functional  run  ex¬ 
plicitly  defining  the  same  k-point  mesh  in  addition  to 
the  desired  k-path  as  specified  in  the  high  symmetry  ex¬ 
amples  of  Ref.  [40] .  The  crucial  step  here  is  to  set  the 
mixing  weights  of  the  extra  k-path  to  zero,  while  keep¬ 
ing  the  original  mesh  intact.  Since  the  orbitals  at  the 
extra  k-path  do  not  contribute  to  the  total  energy,  and 
the  wavefunctions  on  the  original  mesh  are  converged 
as  input,  it  is  only  necessary  to  converge  the  orbitals 
at  the  extra  k-points  mesh  with  the  appropriate  vasp  in¬ 
structions.  The  band  structure  calculation  following  this 
recipe  will  take  approximately  the  same  amount  of  time 
as  a  regular  total  energy  hybrid  calculation  since  the  or¬ 
bitals  at  the  standard  mesh  are  pre-converged. 

As  an  example  Figure  8  shows  the  HSE06  band  struc¬ 
ture  computed  along  the  standard  high  symmetry  fines. 
Comparisons  can  be  drawn  with  the  PBE+U  electronic 
structure  of  Ref.  [87].  The  HSE06  calculation  used 
the  standard  value  of  a  =  0.25.  (Note  that  by  using 
a  =  0.375  the  band-gap  can  be  brought  to  agreement 
with  experiment.)  The  value  of  the  band-gap  is  consid¬ 
erably  improved  Eg  =  2.48  eV  in  comparison  with  GGA 
and  PBE+U  values  of  Eg  =  0.7  eV  and  Eg  =  1.82  eV, 
respectively  (the  phenomenologically  corrected  value 
is  3.36  eV  [86]).  In  other  materials  this  improvement 
is  often  better,  since  ZnO  produces  one  of  the  largest 
LDA/GGA  band-gap  errors.  The  HSE06  band-gaps  for 
narrow  and  medium  gap  semiconductors,  with  the  stan¬ 
dard  ratio  of  the  Fock  exchange,  are  in  good  agree¬ 
ment  with  experiment.  For  wide  band-gap  semiconduc¬ 
tors  and  insulators  the  gaps,  although  underestimated, 
present  a  significant  improvement  over  LDA/GGA  val¬ 
ues  [80].  In  addition,  hybrid  functionals  improve  ef¬ 
fective  mass  estimation  for  almost  all  systems,  typically 
yielding  values  within  a  few  percent  of  experiment  [88]. 

A  fully  functional,  consistent  and  robust  Aflow 
framework  with  hybrid  functionals  is  planned  for  2012. 


7.  Summary 

We  describe  the  Aflow  software  package  for  HT  cal¬ 
culations  of  material  properties.  It  should  be  helpful  to 
the  materials  scientist  seeking  to  determine  properties 


OlZnl_ICSD_26170  (HEX)  Density  of  States  (States/eV) 


Figure  8:  HSE06  hybrid  functional  band  structure  of  ZnO.  The  HSE06 
gap  is  2.48  eV.  LDA/GGA  gap  is  0.7  eV  while  the  PBE+U  value  is 
1.82  eV  [87].  Standard  HSE06  parameters  were  used,  the  ratio  of 
exact  exchange  is  0.25  and  the  screening  parameter  is  0.2  A-1 . 

of  alloy  and  compound  structures.  The  code,  and  op¬ 
eration  manuals  describing  its  features,  are  freely  avail¬ 
able  for  download  at  aflowlib.org/aflow.html.  The  struc¬ 
ture  manipulation  and  analysis  capabilities  of  Aflow 
are  also  available  online  at  the  AcoNVASP-online  entry 
of  this  website. 
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